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Abstract. This paper presents an application of the recent relativistic HLLC 
approximate Riemann solver by Mignone & Bodo to magnetized flows with vanish- 
ing normal component of the magnetic field. The numerical scheme is validated in 
two dimensions by investigating the propagation of axisymmetric jets with toroidal 
magnetic fields. The selected jet models show that the HLLC solver yields sharper 
resolution of contact and shear waves and better convergence properties over the 
traditional HLL approach. 
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1. Introduction 

Jets are one of the fundamental components of radio-loud active galac- 
tic nuclei (AGN). Although different in detail, they all share the prop- 
erty of being relativistic on the parsec scale and of carrying magnetic 
fields that may be dynamically important. Observations of their non- 
thermal synchrotron radiation, in fact, require the presence of a mag- 
netic field and the high degree of polarization observed indicates the 
presence of some large scale structure in the field. 

Theoretical investigations require, therefore, a relativistic MHD de- 
scription of such objects. In this effort, recent numerical simulations 
(Komissarov, 1999b; Leismann et al., 2005) have turned out as ef- 
ficient investigation tools in understanding many aspects of the jet 
evolution. This, in turn, justifies recent and new efforts to extend ex- 
isting gas-dynamical Godunov-type codes to the relativistic regime, see 
Komissarov (1999a), Balsara (2001), Del Zanna et al. (2003), and the 
extensive review by Marti & Miiller (2003). In this perspective, we ex- 
tend the recently derived relativistic hydrodynamical HLLC Riemann 
solver (Mignone & Bodo, 2005) to the magnetic case, with vanishing 
normal component of the magnetic field. The novel numerical scheme 
is then applied to the propagation of relativistic jets in presence of 
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toroidal magnetic fields. The paper is organized as follows: in §2 we 
give the relevant relativistic MHD equations, in §3 we briefly review the 
approximate Riemann solver. Finally, in §4 we discuss the astrophysical 
jet applications. 



2. RMHD Equations 

The motion of a magnetized, ideal relativistic fluid is governed by the 
system of conservation laws (Anile, 1989) 

a 

where d = x,y,z and U = {D, nik, -Bfe, E) is the vector of conservative 
variables with components given, respectively, by the laboratory den- 
sity D, the three components of momentum mk and magnetic field Bj- 
{k = X, y, z) and the total energy density E. are the fiuxes along 
the x'^ direction; for k = x one has 



Dvx, mkVx - Br, + (v • B)v''^ +pSxk, B^Vr, - B^Vk, 

^ (2) 
where p = Pg + |Bp/(27^) + (v • B)^/2 is the total pressure, pg is the 
thermal (gas) pressure and v = {vx,Vy,Vz) is the fiuid three- velocity. 
The Lorentz factor is denoted with 7. Expressions for F^(U) and F^(U) 
follow by cyclic permutations of the indexes. 

Introducing the primitive flow variables V = {p, v, pg,'B), one has 

D = P7, (3) 
ruk = (phj^ + \Bf)vk - (v ■ B)Bk , (4) 

2 |B|2 |v|2|B|2-(vB)2 

where p and h are, respectively, the rest-mass density and specific en- 
thalpy of the fluid. An additional equation of state is necessary to close 
the system (1). Throughout the following we will assume a constant 
T-law, with specific enthalpy given by 

where F is the constant specific heat ratio. 
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2.1. Recovering primitive variables 

Godunov-type codes (Godunov, 1959) are based on a conservative for- 
mulation where laboratory density, momentum, energy and magnetic 
fields are evolved in time. On the other hand, primitive variables, 
V = (p, V, jOg, B), arc required for the computation of the fluxes (eq. 2) 
and more convenient for interpolation purposes. 

This task requires the solution of the nonlinear equation 

fiW)^W-p, + {l-^)\B\^-E-^^=0 (7) 

where W = ph'y'^, S = m • B. Equation (7) follows directly from 
equation (5) together with (4). Since, at the beginning of a time step 
E, B and m are known quantities, both 7 and pg may be expressed in 
terms of W alone from 

|mP = (w.-+|B|f (l-i)-|i(|BP + 2W'), (8) 

and from the equation of state (6) using p = D/^. Equation (7) 
can be solved by any standard root finding algorithm. Once W has 
been found, the Lorentz factor is easily found from (8), the thermal 
pressure from (6) and velocities are found by inverting equation (4): 

Finally, equation (3) is used to determine the proper density p. 



3. A relativist ic HLLC Riemann solver 

In the traditional Godunov approach (Godunov, 1959), the conserva- 
tion laws (1) are advanced in time by solving, at each zone interface 

X - , 1, a Riemann problem with initial condition^ 
'+2 



U(x,0) 



r u 

u 



if X < X, 
if X > X, 



(10) 



i+i 



where U^^ and 1 are the left and right edge values at zone in- 

terfaces. The exact solution to the Riemann problem for the relativistic 



^ In what follows we take the x axis as the normal direction. 
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MHD equations has been recently derived by Giacomazzo & Rezzolla 
(2005). As in the classical case, a seven- wave pattern emerges from the 
decay of an initial discontinuity separating two constant states. The 
resulting structure can be found by iterative techniques and it can be 
quite involved. 

For the present purposes, however, we will be concerned with the 
special case where the component of magnetic field normal to a zone 
interface vanishes. Under this condition, a degeneracy occurs where 
the tangential, Alfven and slow waves all propagate at the speed of 
the fluid and the solution simplifies to a three-wave pattern. In this 
special situation, the relativistic HLLC approximate Riemann solver 
originally derived for the gas-dynamical equations by Toro et al. (1994) 
and recently extended to the relativistic regime by Mignone & Bodo, 
2005 (paper I henceforth), can still be applied with minor modifica- 
tions. This statement stems from the fact that both state and flux 
vectors on either side of the tangential discontinuity (denoted with the 
* superscript) retain the same form as in the non-magnetized case: 



U2,K = [D*, ml, ml, ml, B^, B*,, E*)^^^ , (11a) 

^L,R = {d*vI, mlvl+p*, m*yvl, m*<, "^x)^ ^ , (Hb) 

where we can still set m* = [E* +p*)vl. Total pressure and normal ve- 
locity are continuous across the middle wave, i.e. = p*^, v^. ^ = v^. ^. 
The method of solution is thus entirely anakjgous to the strategy shown 
in paper I and we will not repeat it here. The two additional components 
related to the presence of transverse magnetic field are decoupled from 
the rest and can be solved without additional complications. Once v* 
has been found, in fact, we use the jump conditions to get 

K. = ^^.. (12) 

from the pre-shock states. 

Finally, we need an estimate for the outermost right and left-going 
fast wave speeds, Xr and A^. These two signal velocities may be found 
by exploiting the characteristic decomposition of the RMHD equation, 
extensively analyzed by Anile &; Pennisi (1987) and Anile (1989). In 
the general case, the fast and slow magneto-sonic waves satisfy a quar- 
tic equation solvable by means of numerical and analytical methods, 
e.g. Del Zanna et al. (2003). For vanishing normal component of the 
magnetic field, however, the two roots corresponding to the degenerate 
slow waves have the trivial solution X = Vx, whereas the fast speeds 
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satisfy the following quadratic: 

(A - v:,f = as{l - A^) , (13) 

with 

_ B^/^^ + phcl + {^^-^)\l-cl) 

ph{l - C2)72 ^'^^ 



where = \jTp/ (ph) is the speed of sound. The solutions of (13) are 
used to provide suitable guesses for the outermost signal velocities in 
our HLLC Riemann solver. As in paper I, we set 

Xl = min(A_(VL),A_(Vj?)) , Xr = max (A+(Vl), A+(Vfl)) , (15) 

where A_ (A+) is the smallest (biggest) fast wave speed. 



4. Code Validation 

Spatial and temporal discretization of equations (1) is based on the 
second-order Hancock predictor step already described in paper I. For 
the sake of conciseness we will not repeat it here. 

4.1. A Shock-Tube Example 

As an illustrative test, we consider an initial discontinuity separating 
two constant left and right states given by 

r (1, 30, 0, 0, 0, 20, 0)^ for x < 0.5 , 

(16) 

The domain is the interval [0, 1] covered with 1600 uniform compu- 
tational cells. The ideal equation of state (6) is used with T = 4/3. 
Figure (1) compares the result computed with the HLLC solver with the 
analytical solution (available from Giacomazzo & RezzoUa (2005)), 
at the final time t = 0.25. The breakup of the discontinuity results 
in a left-going fast rarefaction wave and a thin shell of high density 
material bounded by a tangential discontinuity and a gasdynamical 
shock propagating to the right. The total pressure and magnetic field 
are continuous through the contact and the shock wave, respectively. 
All discontinuities are correctly captured, with satisfactory resolutions: 
the tangential discontinuity and the shock smear out over ~ 5 -^ 6 and 
4-^5 zones, respectively. Small overshoots appear at the tail of the 
rarefaction fan. A similar, perhaps more pronounced behavior is shown 
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Figure 1. Computed profiles (symbols) and exact solution (solid line) for the 
first test problem. The second-order scheme with the HLLC Rieman solver has 
been employed on 1600 zones. The final integration time is t = 0.25; the Courant 
number is 0.5. The solution is comprised of a loft-going fast rarefaction, a right-going 
tangential and shock waves moving close to each other. 

in the results by Komissarov (1999a) who used a linearized, Roe-like 
Ricmann solver. The relative error in the thin shell is < 1% for both 
pressure and density. 



4.2. AXISYMMETRIC JET WITH TOROIDAL MAGNETIC FIELD 

We now turn our attention to the propagation of a relativistic magne- 
tized jet in cylindrical coordinates {r,z). For the sake of comparison, 
we have implemented three models already discussed in literature. 
The first setup, denoted with 5*1, is model A of Komissarov (1999b) 
(K99 henceforth), the second and third setups, 52 and S3, are taken 
from models C2_10/3 and C2_l/20 of Leismann et al. (2005) (L05 
henceforth), respectively. 

In all models, the static external medium has uniform density Pe and 
zero magnetic field. At the inlet, r < 1 and z = 0, a supersonic beam 
with uniform velocity Vj (Lorentz factor 7^) and density pj is injected 



rmhd.tex; 5/02/2008; 9:00; p. 6 



RMHD Simulations 



7 



Table I. Relevant parameters for models SI, 5*2 and S3 de- 
scribed in the text. Prom left to right: jet density, ambient 
medium density, beam velocity, magnetization radius, average 
magnetization, specific heat ratio and magnetic to rest mass 
energy ratio. 



Model 


Pi 


Pe 






P 


r 


a 




SI 


1 


10^ 




0.37 


1 


4/3 


0.34 


S2 


10-2 


1 


0.99 


0.6 


10/3 


5/3 


0.11 


S3 


10-2 


1 


0.99 


0.6 


1/20 


5/3 


0.0017 



with a toroidal magnetic field obeying 

r/vm if r < 
Ijhmrm/r if r^^ < r < 1 , 



Mr) 



(17) 



where is the magnetization radius of the beam. The ideal equation 
of state (6) is used in all calculations. In model 51, we set 6^ = 1 
and the thermal pressure inside the beam, pj, is prescribed from the 
condition of hydromagnetic equilibrium, 



Pj{r) = < 



Pe 



a + I (1 - {r/r„ 



ape 



if 

if 



r <rm, 
rm <r <1 , 



(18) 



where a = 1 — r^/ f3, /3 = 0.34 and pe = /S^m/^ is the thermal pressure 
of the ambient mediTim. In models 5*2 and S3 bm is computed from the 
average magnetization parameter (3: 



rl^{l - 41ogr„ 



(19) 



and the thermal pressure follows directly from the definition of the 
Mach number, i.e. 



Pj 



Pi^l(r-i) 



r(r - i)M/ - Tv' 



(20) 



with Mj = 6. The numerical values of the relevant parameters are 
given in Table I. The simulations are performed on the physical domain 
r G [0, 10.5], z G [0,50] using 20 points per jet radius. The total grid 
size is 210 x 1000. 
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Figure 2. Model SI at t = 100. The top and bottom panels show, respectively, the 
density logarithms and magnetization parameter (5 for the HLLC (upper half of 
each panel) and HLL (lower half of each panel) runs. Notice that the two panels use 
opposite color gradients. The grid resolution is 20 zones per jet beam, giving a total 
grid size of 210 (in r) and 1000 zones (in z). Integration has been carried with CFL 
= 0.4 and the second-order Hancock scheme outlined in paper L 



For each model, we have carried two set of simulations: the first one 
employing the HLLC approximate Riemann solver described in §3 and 
the second one adopting the HLL scheme, also used in L05. 

Models 51 and S2 are characterized by a high Poynting flux B'^v 
and tend to develop prominent nose cones due to the strong magnetic 
pinching. This is evident from Fig. (2) and (3), where the upper (lower) 
portion of each gray scale image refers to the HLLC (HLL) approximate 
Riemann solver at t = 100 and t = 126, respectively. The nose cone 
appears as the extended high pressure region bounded by the terminal 
conical shock (mach disk) and the bow shock. At the Mach disk, the 
beam is strongly decelerated and magnetic tension wards off sideways 
deflection of shocked material into the cocoon. The magnetization pa- 
rameter, defined as /3 = i?|/(27^pg), is particularly strong in the cone 
above the symmetry axis. As pointed out by K99, morphological prop- 
erties are prevalently dictated by the ratio of the kinetic energy flux to 
the Poynting flux, k = pjhjjj / B^, which determines the magnetization 
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Figure 3. Model S2 at t = 126. The top panel shows the density logarithm in 
the HLLC (upper half) and HLL (lower half) runs. The bottom panel shows the 
magnetization parameter in the HLLC (upper half) and HLL (lower half) runs. The 
final integration time is the same one as in L05. 



of shocked gas behind the terminal shock. An alternative parameter 
(L05) is the ratio of the magnetic energy density to the rest mass, 
a = i?^/(7^p): extended nose cones tend to form for low k (high a), 
whereas familiar cocoon structures appear for high values of k (low 
a), regardless of the jet magnetization parameter f3. Values of a are 
reported in Table I. The positions of the reverse shocks for the HLL 
and HLLC runs are, respectively, z ^ 30 and z 26 in model S*!, z ~ 33 
and z ~ 41 in model S'2. However, the most striking difference between 
the HLL and HLLC schemes is the jet propagation speed obtained from 
model SI. Indeed, at t = 100 the head position in the HLLC run is 
z ~ 47.6 to be compared with a value of z ^ 41.8 obtained from the 
HLL integration. Note that the final integration time is not the same 
as in K99, where the simulation ends at a later time, t = 110, yielding 
a head position value of z 46, apparently more consistent with the 
HLL run. However, in order to further investigate these discrepancies, 
we carried out an additional set of simulations at half and twice the 
resolution. Figure (4) plots the jet head position as function of time 
for the three resolutions in consideration: low (105 x 500), medium 
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50 F 




time 

Figure 4- Jet head positions for model 51 computed for three different resolution 
runs: low (105 x 500), medium (210 x 1000) and high (240 x 2000). The HLLC runs 
are shown using, respectively, dotted, dashed and solid lines; for the HLL we use 
plus signs, triangles and squares. 

(210 X 1000) and high (420 x 2000). The HLL scheme achieves the same 
convergence of the HLLC scheme only at high resolution, whereas the 
propagation speeds in the HLLC runs show to be by far less sensitive 
to the grid size. 

This disagreement is remarkably reduced in the second model 5*2, 
where the HLL runs shows a slightly higher propagation speed, but 
again more evident in model 5*3, although with an opposite trend. Once 
more we ascribe this behavior to resolution effects, based on the fact 
that the same model in L05 with twice the resolution (i.e. 40 zones 
per beam radius instead of 20) shows stronger affinities with the HLLC 
integration presented here. 

When the HLLC solver is employed, the backflow turbulent patterns 
in the cocoon are better defined than with the HLL scheme, where the 
larger numerical viscosity inhibits the formation of small scale structure 
at lower resolution. This is corroborated by a simple qualitative com- 
parison with the results of LOS, who used twice the resolution employed 
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Figure 5. Model S3 at t = 126. The final integration time is the same one used in 
LOS. The meaning of the gray scale images in the same as in Fig. (2) and (3). 



here and the HLL solver. The small scale structural details provided 
by the HLLC approach appear to be almost equivalent even at half the 
resolution. This consideration together with the convergence test pre- 
viously shown confirms, indeed, that the ability to describe tangential 
and shear waves largely improves the quality of solution, particularly 
in the multidimensional case. 

Finally, we remind that the computational cost of the HLLC scheme 
over the traditional HLL Ricmann solver is less than 7% and we believe 
that this largely justifies its use in Godunov-type codes. 



5. Discussion 

Numerical simulations of relativistic magnetized jets with toroidal mag- 
netic fields have been presented. The three models considered show 
formation of conspicuous nose cones and jet confinement with increas- 
ing magnetization. The morphological properties favorably compare 
with the results of previous investigators. Simulations have been con- 
ducted by using the recent relativistic hydrodynamics HLLC solver by 
Mignone k, Bodo, opportunely extended to the special case of vanishing 
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normal component of magnetic field. In order to assess the validity of 
the new algorithm, simulations have been directly compared with the 
traditional HLL Ricmann solver (Del Zanna et al., 2003). It is found 
that the HLLC approach largely improves the quality of solution over 
the HLL scheme by considerably reducing the amount of numerical vis- 
cosity. This is manifested in a sharper resolution of fine scale structures 
in proximity of shear and tangential flows. Furthermore, convergence 
tests at different resolutions have shown that, with the HLLC scheme, 
the jet head position has little dependence on the grid size. This is 
a remarkable property for a shock-capturing scheme and has to be 
consolidated through further numerical testing. Generalization to the 
full relativistic MHD case is underway. 
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